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On the physical origin of dark matter density profiles 
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ABSTRACT 

The radial mass distribution of dark matter haloes is investigated within the frame- 
work of the spherical infall model. We present a new formulation of spherical collapse 
including non-radial motions, and compare the analytical profiles with a set of high- 
resolution N-body simulations ranging from galactic to cluster scales. We argue that 
the dark matter density profile is entirely determined by the initial conditions, which 
are described by only two parameters: the height of the primordial peak and the 
smoothing scale. These are physically meaningful quantities in our model, related to 
the mass and formation time of the halo. Angular momentum is dominated by velocity 
dispersion, and it is responsible for the shape of the density profile near the centre. 
The phase-space density of our simulated haloes is well described by a power-law pro- 
file, p/er 3 = 10 1 - 46±0 04 (p c /w 3 ir )(r/r vir )- 1 - 90±0 05 . Setting the eccentricity of particle 
orbits according to the numerical results, our model is able to reproduce the mass 
distribution of individual haloes. 
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1 INTRODUCTION 

The hierarchical clustering paradigm states that the growth 
of cold dark matter (CDM) haloes proceeds by accretion of 
smaller units from the surrounding environment, either by 
continuous infall or by discrete merging events. 

Cosmological N-body simulations are a valuable tool 
to study the mass distribution of dark matter haloes and 
its evolution in the non- li n ear re gime. In the early numer - 
ical work of iQuinn et alJ il986T) and iFrenk et alJ lll98St) , 
haloes showed an is o therm al de nsity profile ( p oc r~ 2 ). 
iDubinski fc Carlberd dl99lT) and ICrone et all <199-fl had 
enough resolution in their simulations to detect the first 
evidence of d e parture fro m a pure power-law. Later on, 
iNavarro et"al] Il996t Il997t hereafter NFW) found that the 
density profile could be fitted by a simple analytical function 



scales could be described in terms of a one-parameter family 
of analytical profiles. 

Similar results have been found in independent simula- 
tions with much higher mass and force resolution than the 
original NFW paper. However, there is still some controversy 
about the innermos t value of a a nd its dependence on reso- 
luti on. iMoore et alJ J199& ll99Sft . iGhigna et al.1 Jl998t l200d) 
and iFukushige fc Making (119971 . 12001) find steeper density 
profiles near the centre (a ~ —1.5), w hereas other authors 
bine fc Sutoi2000l:lKlvpin et all200j) claim that the actual 
value of a may depend on halo mass, merger history, and 



value or a may depend or 
substructure. iPower et alJ 



p(r) 



(r/r s )(l + r/r s 



(1) 



in terms of a characteristic density p B and a characteristic 
radius r s . This profile is steeper than isothermal at large 
radii and shallower near the centre. The logarithmic slope 
of the density profile, a(r) = dlog(p)/dlog(r), tends to a = 
—3 for r — ■» oo and a = — 1 for r — > 0. It corresponds to 
the i sothermal case , a = —2, at the characteristic radius 
only. INavarro et alJ il997fl further showed that the two free 
parameters in equation Q are not independent. Should this 
be true, the final mass distribution of objects of different 



2003) pointed out that the log- 
arithmic slope becomes increasingly shallow inwards, with 
little sign of approaching an asymptotic value at the re- 
solved radii. In that case, the precise value of a at a given 
cut-off scale would not be parti cularly m eaningf ul. Th is re- 
sult has been later co nfirmed bv IFukushige et al.1 J2003I) and 
lHavashi et alJ <|2003T). and it is predicted by several analyt i- 
cal models fe.g. lTavlor fc NavarrolEoOil iHoeft et alJl2003ft . 

Observed r otation curves of d warf s piral an d LSB 
galaxies (e.g. iFlores fc Primackl 11994 iMoorel Il994t 
[ Burked Il995l: iKravtsoy et alJ I l998fc iBorriello fc Saluccil 
20011 Ide Blok et al.l l200lt Ide Blok fc Bosmal 120021 
iMarchesini et ahl I2002T) seem to indicate that the shape 
of the density profile at small scales is significantly 
shallower than what is found in numerical simulations. 
This discrepancy has been often signalled as a genuine 
crisis of the CDM scenario, and several alternatives 
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have been suggested, su c h as warm jColm et al-j |2000l: 
I Sommer-Larsen fc Dolgovl l200ll) . repulsive l|Goodmanl 
l200(t>. fluid ]Peeblesll200rl . fu zzv JHu et alJl20Qnh. de cay- 
ing llCenl 1200 IT ), ann ihilating (Kaplinghat ct al. 200 0J), or 
self-interacting llSpergel fc Steinhlu^ltil2'o'o'oT ~ [Yo"sIiida et alJ 
" iDave et aljl200lft dark matter. 



Unfortunately, it has proved remarkably hard to es- 
tablish the inner slope of the dark matter distribution 
observ ationally (see e. g. |SwMers_et_aj |2003l) . Some au- 
thors l|v an den Bosch fc Swater^l2^0^njimenez et ai1l2003t 
ISwaters^tHnTl^OO.?!) claim that a cuspy density profile with 
a ^ — 1 is consistent with current observations, although a 
shallower slope is able to explain them as well. Yet, a value 
as steep as a = — 1.5 can be confidently ruled out in most 
cases. According to lHavashi et alJ ]2003l) . only about 30 per 
cent of the rotational curves are actually inconsistent with 
the simulation data. 

On cluster scales, X-ra y analyses have led to wide rang- 
ing results, f r om a = -0.6 fcttori et al.ll2002ft to a = -1.2 
llLewis et al] l2003l) or even a = —1.9 jArabadiis et al] 
l2002T) ~Measurements based on gravitational lensing yield 
conflicting estimates as well, either in rough agreement 
with the results of nume rical simulations (e.g. Dah le et al] 
l2003l : lGavazzi et al.||2003 |). or finding m uch shallower slopes, 
a = -0.5 (e.g. lSand et al]l20oll2003h . 

A conclusive theoretical prediction of the central mass 
distribution of CDM haloes is therefore an important check 
for any model of structure formation. The controversy re- 
garding the 'universal' density profile and its logarithmic 
slope at the centre has stimulated a great deal of analytical 
work. On one hand, we would like to find out not only the 
actual shape of the profile, but also the physical mechanisms 
behind the 'universality' observed in numerical N-body sim- 
ulations. On the other hand, it would be interesting to find 
and explain correlations with halo size, environment, power 
spectrum or even the nature of dark matter particles. 

A number of plausible arguments about the radial struc- 
ture of CDM haloes have been advanced during the last 
30 years. The basic problem of the collisionless collapse of 
a spherical perturbation in an expanding b ackground was 
first a ddre ssed in the tw o seminal papers by iGunn fc Gott 
]l972l) and lGunnl ]l977l) . where the cosmological expansion 
and the role of adiabatic invariance were first introduced 
in the context of the formati on of individual objects. The 
nex t step was acc o mplis hed bv lFillmore fc Goldreichl (I1984T) 
and iBertschingerj 1^)85), who found analytical predictions 
for the density of collapsed objects see ded by scale-free pri- 
mordi al perturbations in a flat universe. iHoffman fc Shahaml 
(1985) generalised these solutions to realistic initial con- 
ditions in flat as well as open Friedmann models. Modifi- 
cations of the self-similar collapse model to include more 
realistic dy namics of the growth p r ocess have been pro- 
posed (e.g. lAv ila-Reese et al.| |l998fr |Henriksen fc Wid rowl 
ll999l:lLokasi200d : lKulll999llSubramanian et all2000l). Sev- 
eral a uthors (e.g. |Sj[e rfcWhite] [1998|; [Salvador- S ol e et al] 
119981: INusser fc Shethl Il999l : iManriaue et all l2003Tl argue 
that the central density profile is linked to the merging his- 
tory of dark matter substru cture, and baryons have been 
invoked bo th to shallow (e.g. | El-Zant et al.ll200ll 120031) and 
to steepen ]Blumenthal et ah ^86Ti ~tIie dark matter profile. 

In this paper, we present an analytical model for 
the assembly of CDM haloes based on the spherical col- 



lapse paradigm. Following the spirit of IHoffman fc Shahaml 
( 1985), we will assume that objects do not form around local 
maxima of the primordial density field, but of the smoothed 
density field. We argue that all information about the ini- 
tial conditions below the smoothing scale is lost during the 
merging process. We will show that setting the smoothing 
scale for a halo of a given mass is equivalent to specify- 
ing its formation time. Initial conditions are computed fol- 
lowing the Gaussian random peaks statistics described by 
iBardeen et al] dl986l hereafter BBKS), and angular momen- 
tum is included in a phenomenological way. Our theoreti- 
cal predictions are then compared with the results of high- 
resolution numerical simulations, showing that this model is 
able to accurately reproduce the mass distribution of indi- 
vidual objects. 

The paper is structured as follows. Details of our imple- 
mentation of spherical collapse are given in Section [5] Nu- 
merical experiments are described in Section |H] where both 
the mass and velocity distributions are investigated. We dis- 
cuss our results in Section 2] and Section summarises our 
main conclusions. 



2 SPHERICAL COLLAPSE 

The assembly of dark matter haloes is a highly non-linear 
process, and strong simplifying assumptions must be made 
in order to tackle the problem analytically. Traditionally, 
there a re two complementar y paradigms: the spherical infall 
model jGunn fc GottHl972l)' and the Press-Schecter formal- 
ism JPress fc Schechteill97l). in which mergers play a dom- 
inant role (see e.g. INusser fc Shethl Il999t IManriaue et al] 
2003). 



2.1 The model 

The most simple way of addressing the problem of structure 
formation is t o ass ume spherical symmetry. As shown by 
iTolmanl Jl934 and iBondil lll947l) . a spherically symmetric 
solution of the Einstein equations can be easily interpreted 
in terms of Newtonian dynamics. The equation of motion 
for a Lagrangian shell enclosing a mass M can be derived 
from the conservation of energy 

m 2 r 
where r(t) is the Lagrangian coordinate, M(r) — Mi(n) 
is the mass contained within the shell, and pA = is 
the vacuum energy density associated to a cosmological con- 
stant. Throughout this paper, the subscript 'i' will be used 
to denote initial conditions. 

If we assume the universe to be homogeneous, the mass 

is given by Mi — ^f2J n Pc r ? , where p' c = Substituting 
in @ and making a coordinate transformation 

r(t) = na(t) (3) 

we obtain the Friedmann equation without radiation: 



H? 



fl\a = 1 — Q % „ 



(J) 



In a homogeneous universe, a(t) is independent of n 
and plays the role of a uniform cosmic expansion factor. 
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Our model assumes that structures grow from spherically 
symmetric perturbations, defined as 

M(n) 



Ai(n) 



- 1 



(•>) 



with < A; <C 1. These perturbations are introduced at an 
early epoch ai by slightly displacing the spherical shells of 
matter. Keeping only terms to first order in Ai(n), the new 
positions and velocities are given by 



Hm ( 



1 - -Ai(n 
3 u 



(6) 



The mass enclosed by the perturbed shell is still Mi, 
and the initial specific energy (using ~ 1 and Q,\ ~ 0) 
is £i ~ — |Ai(//jri) 2 . Energy conservation © leads to the 
equation of motion 
• 2 

(7) 



5 . a~ 



According to this equation, the evolution of a single 
spherical shell would be similar to that of a Friedmann uni- 
verse. For high enough values of Ai, the shell reaches a 
maximum radius r m at a turn-around time t m and then 
re-collapses. In an Einstein-de Sitter universe (fi m = 1, 
SI a = 0), equation J7J can be solved analytically, yielding 



3n 
5Ai 



(8) 



2f/i(5Ai/3) 3 / 2 ' 

This is also a valid approximation for shells reaching 
their turn-around radius r m before the cosmological constant 
term starts to dominate the expansion. Since the shells need 
at least another i m to virialise, expression JSJ can be applied 
in a ACDM universe to estimate the maximum expansion 
radius and time for the innermost part of a virialised halo. 
For the outer shells, the equation of motion J7J must be 
integrated numerically in order to find the trajectories up 
to the maximum radius. 

In the absence of shell-crossing, shells would reach the 
origin at T = 2t m . Since they are assumed to be composed of 
collisionless CDM particles, they would simply pass through 
the centre, describing an oscillatory motion with amplitude 
r m and period T. However, equation Q holds as long as the 
enclosed mass Mi remains constant. As a shell re-collapses, 
its particles will cross the orbits of inner shells, and energy 
will not be a constant of motion any more. 

After turn-around, our model assumes that the parti- 
cles belonging to a shell oscillate (or, more generally, orbit) 
within the gravitational potential of the dark matter halo. 
Since CDM particles are expected to spend most of the time 
in the outermost regions of their orbits (particularly when 
angular momentum is taken into account), we approximate 
the mass distribution of the halo by a simple power law 

M(r) = Mi ( — (9) 
V r m / 

where qmH = d log M(r)/d log r is the local value of the 
logarithmic slope of the mass profile, evaluated at the max- 
imum radius of the orbit. 

At first sight, this might see m similar to the cla ssical 
approach based on self-similarity feertschingerl ll985). but 
in that case the final mass distribution is indeed assumed 
to be a power law, whereas in our model this ansatz is only 
an approximation to compute the local potential. The final 



mass profile is obtained self-consistently, adding the contri- 
butions from all shells, each of them with an individual value 
of o M (r m ). 

The probability of finding a particle inside radius r is 
proportional to the fraction of time it spends within r: 



P(r, r n 



da; 
v r (x) 



1 



dx 



^<E>(r m ) - <S>(x) 



(10) 



We evaluate numerically the value of OM(r m ) in order to 
compute the potential. Taking different prescriptions, such 
as an isothermal profile (ckm = 1 for every shell) or a Kep- 
lerian potential (ckm = 0) does not lead to qualitative vari- 
ations in the probability P(r, r m ) and the resulting mass 
distribution. 

If phase mixing is considered to be efficient, particles 
initially on the same shell will be spread out from r — 
to r = r m a short time after i m . For the sake of simplicity, 
we will consider that phase mixing is instantaneous, so im- 
mediately after turn-around the shell is transformed into a 
density distribution whose cumulative mass is proportional 
to P(r,r m ). 

This means that recently accreted particles contribute 
to the mass within x m < r m (i.e. shell-crossing). For shells 
whose maximum radius was x m , the enclosed mass changes 
from Mi(xm) to 



M(x n 



Mi(Xxa) +M add ( 



(11) 



where M a dd(2?m) accounts for particle s belonging to outer 
shells . To compute M a dd(a;m) (see IZaroubi fc Hoffman! 
1993), we must integrate the contribution of every shell 
whose maximum radius is larger than x m , up to the cur- 
rent turn-around radius i? ta : 



Madden 



ta dMi(r) 
dr 



P(x m ,r) dr. 



(12) 



To compute the evolution of th e shell afte r shell- 
crossing, we apply adiabatic invariance jGunnlll977l) . If the 
potential evolves on a timescale much longer than the or- 
bital period of the inner particles, their dynamics admits an 
adiabatic invariant 



Jr = i> v r( r ) dr 



(13) 



also known as the radial action. For a power-law potential, 
the radial action J r is proportional to y x ul M{x ul ). When 
we increase the mass by an amount M a dd, the maximum 
radius of the inner shell must decrease in order to keep J r 
constant. The final radius, xo, is given by the implicit equa- 
tion xo = F(xo) Km, where 

Mi Mi 



F(x ) 



(14) 



Mi + M add (xo) Mi + M add [F{x )x m ] 

and whose solution must be obtained numerically for each 
shell. 

To summarise, the numerical procedure to compute the 
final radius ro(to) of a Lagrangian shell of matter involves 
the following steps: 

(i) Set Mi (or, equivalently, n). Start by the outer shell. 

(ii) Integrate the equation of motion J7J up to t m . 

(iii) If t m > to, the shell is still expanding: ro = r(to). 

(iv) If t m < to, solve 1141 to compute ro = -F 1 (ro)r m , and 
add the contribution of this shell to M a dd(f). 
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(v) Repeat for the next shell towards the centre. 



2.2 Initial conditions 

The model described above allows us to compute the mass 
distribution arising from a primordial fluctuation Ai(n), but 
does not say anything about the shape of this function or 
its physical origin. Never the less, it is important to note 
that the final density profile is entirely determined by this 
initial condition. Thus, in the spherical collapse paradigm, 
the case for a universal density profile can be reformulated 
in terms of universality in the primordial fluctuations that 

set Ajfn). 

iHoffman fe Shahaml lll985f) suggested that, according 
to the hierarchical scenario of structure formation, haloes 
should collapse around maxima of the smoothed density field. 
The statistics of peaks in a Gaussian random field has been 
extensively studied in the classical paper by BBKS. A well- 
known result is the expression for the radial density profile 
of a fluctuation centered on a primordial peak of arbitrary 
height: 



(*(r)> 



O~0 



7(1 -7 2 ) 



7 2 T/i(r) + ^V 2 ^(r) 



(15) 



where i/>(r) = £(r)/ao is the normalised two-point corre- 
lation function, <to = £(0) 1,/2 is the rms density fluctu- 
ation, and i/<7o is the height of the peak. The quantities 
7 = t7 2 /(<T2(To) and R, = \/3o"i/o"2 are related to the mo- 
ments of the power spectrum, 



dk, 



(16) 



and the function #(7, 7^) parametrizes the second derivative 
of the density fluctuation near the peak. BBKS suggest the 
approximate fitting formula 

,(7,7.)^ 3(1 : 72) + (L216 - 0V) ^ him2] (17) 



[3(1 



valid to 1 per cent accuracy in the range 0.4 < 7 < 0.7 and 
1 < 7^ < 3, which is the scale relevant for both galaxies and 
galaxy clusters. 

Expression <!15t is often quoted in t he literature (e.g. 
Hoffman! Il988l: lLokas fc Hoffman! l2000t iDel PopoIo et alJ 
20001: lHioteli3l2002f) as the initial condition Aj(n). However, 
we argue that (<5f(r)} denotes the smoothed density profile 
around a local maximum of the smoothed density field, 



(18) 



Sc(r) = / Wi(r~x)S(x)dx 



where the function Wiif — x) is a smoothing kernel that 
depends on a certain filtering scale Rf. 

The smoothed profile (5f(r)) given by 1151 is in general 
not equal to the mean value of the actual overdensity, <5(r), 
that must be integrated to compute Aj(n): 



Ai(n) = 4-7T / 8{x)x Ax 



(19) 



Comparing this expression with 11811 . we see that Ai(n) 
is equivalent to <5f (0) as long as W{ is taken to be a spherical 
top hat of radius rj. 



To sum up, we are interested in the physical density pro- 
file Ai(n) around a local maximum of the smoothed density 
field. To locate the maximum, we use a Gaussian smoothing 
kernel 



W f (r- 



{2ttR 2 



-3/2 



exp 



2R 2 



(20) 



in order to avoid the oscillations in Fourier space that arise 
from a top hat filter. We set the scale of the fluctuation, R{ , 
and impose the condition that r = is a maximum of 5f(f). 

Then, we compute Ai(n) = S Ti (0) by applying a top 
hat smoothing of radius n. BBKS show that the probability 
distribution of 8 Ti (0) is a Gaussian with mean 



(5 ri (0)) = Vi- 
vo f 

The moments 



1 - 7 2 0-0/ 



'lh "Of 
J 0h °tf 



2(i+i) 



are computed from 
P f (fc) = P{k) exp [-(kR { ) 
and 

P h (k) = P(k) exp {kRc)2 



dk 



(21) 



(22) 



(23) 



3[sin(fcri) — fen cos(fcn)] 
(fen) 3 



(24) 



where P(k) is the ACDM power spectrum that we used to 
generate the initial conditions for our simulations, evaluated 
at time ai. 



2.3 Angular momentum 

Two d ecades after the semina l paper by iGunn fc Gotd 
Jl972h . lwhite fc Zaritskvl Jl992l) pointed ou t (see also the 
pioneering work of iGurevich fc ZvbirJ H988T) that angular 
momentum would prevent the orbits of CDM particles from 
reaching the origin. 

For pure radial orbits, the mass in the centre is dom- 
inated by M a( jd (i.e. particles from the outer shells) when 
the profile at turn-around is shallower than isothermal 
jFillmore fc Goldreichlll984h . The final density distribution 
is p{r) oc r~ 2 , independent o n the initial logarithmic slope . 
More recently, lLokasI fcOOOl) and lLokas fc Hoffmanl J200(J) 
found a similar result considering non-self-similar primor- 
dial fluctuations based on the peak formalism. 

Angular momentum arises in linear theory dWhiteil984h 
from the misalignment between the asymmetry of the col- 
lapsing region (i.e. the inertia tensor) and the tidal forces 
it experien ces during the linea r regime. In addition, violent 
relaxation jLvnden-Belllll967D transforms the radial infall 
energy into random velocity dispersion, which has a tangen- 
tial component. 

The total amount of angular momentum acquired by 
the sys tem is, however, an open question. The usual ap - 
proach jAvila-Reese et al.lll998l : lNusserl200ltlHiotelis)l2002l) 
consists in assigning a specific angular momentum at turn- 
around 



j oc VGMr m . (25) 
With this prescription, the orbit eccentricity e is the 
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Figure 1. Radial density profile arising from a 3cr fluctuation on 
1 h~ 1 Mpc scales (solid lines), for different values of the orbit 
eccentricity e. A NFW profile (dotted line) is plotted for compar- 
ison. 



same for all particles in the halo (lNusserll200lT) . The peri- 
centric radius is given by 

1 — e 

fmin = ~ : ?"max (26) 

1 + e 

where r max is computed from adiabatic invariance 1141 . fol- 
lowing the procedure explained in Section \l. II Angular mo- 
mentum is taken into account by adding the usual term 
j 2 /(2r 2 ) to the gravitational potential $(r) in equation HUH . 
Although this changes the actu al value o f the radial ac- 
tion, J r is still proportional to ^/r m A4"(r m ) and l|14[l can be 
used to compute the collapse factor F(r). The assumption 
of spherical symmetry leads to angular momentum conser- 
vation, which implies constant e during the contraction. 

The effect of angular momentum is illustrated in Fig- 
ureQ We plot the mass distribution at z = corresponding 
to a 3cr peak in the primordi al density field, smoothe d on 
1 ft" 1 Mpc sca les. As noted bv lWhite fc Zaritskvl il992f) and 
iHiotelis) i2002l) . angular momentum plays a key role during 
secondary infall, decreasing the amount of mass contributed 
by recently accreted shells in the innermost regions of the 
halo. 

In agreement with lHiotelij i2002l) . we find that the pre- 
dicted density profile becomes shallower in the centre as the 
amount of angular momentum is increased (lower e). Pure 
radial orbits give rise to a steep profile, so mewhat similar to 
the form proposed bv lMoore et alJ <ll999l) . although the ex- 
act shape is slightly different. When the eccentricity is low, 
the halo develops a constant density core inconsistent with 
the results of numerical simulations. A value e = 0.5 yields a 
mass distribution virtually indistinguishable from the NFW 
formula (dotted line in Figure down to kpc scales. How- 
ever, the logarithmic slope predicted by the model keeps in- 
creasing towards the centre, while NFW tends to an asymp- 
totic value, a — — 1. 



3 NUMERICAL EXPERIMENTS 

In order to test the analytical model, we carried out a se- 
ries of high-resolution N-body simulations of cluster for- 
matio n with the adaptive mesh code ART dKravtsov et alJ 
119971) . Slightly lower spatial resolution experiments have 
been run from the same initial cond itions with the Tree- 
SPH gasdynamical code Gadget dSprineel et alJ 1200 ll : 
ISprineel fc HernauisdEool) . The radial structure of both 
dark and bary o nic c omponents has been addressed in 
lAscasibar et al.l i2003t) . For a detailed descript ion of the 
numer ical experiments, the reader is referred to lAscasibarl 
(2003). 

A sample of cluster-size haloes has been selected from 
an initial low-resolution (128 3 particles) pure N-body simu- 
lation of a 80 hT 1 Mpc box in a ACDM universe (f2 m = 0.3; 
Qa = 0.7; h — 0.7; as = 0.9). Higher resolution has 
been achieved by m eans of the multiple-mass technique (see 
iKlypin et aT]l200lj . for details), using 3 levels of mass refine- 
ment. This is equivalent to an effective resolution of 512 3 
CDM particles (3.1 x 10 s h" 1 M ) in the highest refine- 
ment level. The minimum cell size allowed in the ART runs 
was 1.2 ft -1 kpc. All simulations were started at z = 50. 

This procedure has been applied to 15 objects, ranging 
from 3 x 10 13 to 2 x 10 14 h' 1 M Q . In order to explore a 
broader mass range, a smaller box (25 Mpc) has been 
simulated, and six galaxy-size haloes have been added to the 
cluster sample. Mass resolution is 1.2 x 10 6 ft -1 Mq for these 
objects, with a minimum cell size of 0.2 h~ x kpc. 

All CDM haloes have been classified according to 
their dynamical state according to a substructure-based cri- 
terium. We use th e Bound Density Maxima galaxy fin ding 
algorithm (see e.g. lColfn et aljll999] : lKlvpin et alJll999l) and 
look for the most massive subhalo within the virial radius. 
We label as major merger any object in which the subhalo 
is heavier than 50 per cent the mass of the main object; if 
the mass is between 10 and 50 per cent, it is classified as 
a minor merger; otherwise, we assume the object to be a 
relaxed system in virial equilibrium. 

3.1 Mass distribution 

The spherically-averaged mass distribution of our numerical 
CDM haloes is shown in Figure All profiles have been 
rescaled by their best-fitting characteristic density and ra- 
dius. As can be seen in the Figure, the NFW formula pro- 
vides a reasonable approximation to the mass distribution, 
although some deviations occur at the innermost part, as 
well as beyond the virial radius. 

In Figure [3J the logarithmic slopes of both mass and 
density profiles are plotted according to our dynamical clas- 
sification. We do not find any evidence of an asymptotic 
behavior up to the reso lution limit, in agreement with re- 
cent numerical studies (|Power et al ] l2003l : iFukushiee et alJ 
120031: lHavashi et alJl2003l) . An extrapolation of our results 
suggests that the central slopes in relaxed haloes could 
be indeed less steep than the NFW fit, as predicted by 
analytical models based on the velocity dispersion profile 
jTavlor fc NavarrcfeOOltlHoeft et aT]|2003l) . More resolution 
is clearly needed in order to establish a firm c onclusion. 

As pointed out by s everal authors (e.g. | Jhig fc Sutol 
120001: IKlypin et all 1200 it lAscasibar et alJl2003l) . the mass 
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Figure 2. Radial density (left) and mass (right) profiles of our 
dark matter haloes, rescaled by their best-fitting NFW parame- 
ters and plotted as dotted lines. The corresponding logarithmic 
slopes are shown on the bottom panels. NFW model is indicated 
by a solid line. 



Figure 3. Logarithmic slopes of the mass (top) and density (bot- 
tom) profiles. Black squares show the average over relaxed haloes. 
Minor and major mergers are represented by empty squares and 
stars, respectively. Solid lines indicate NFW model. 



distribution near the centre might depend on the dynamical 
state of the halo. At the resolution of the present simula- 
tions, relaxed haloes are well described by NFW model, but 
merging systems display steeper profiles. A conjecture that 
deserves further investigation is whether violent relaxation 
and subsequent equipartition of energy may drive the in- 
ner part of the system into a nearly isothermal state. That 
would provide a natural mechanism for the final profile to 
be independent on the details of the merging history. 

Another difference between merging and relaxed sys- 
tems is that mergers are typically more extended than re- 
laxed haloes of the same mass. We plot in Figure 0] the 
mass-con centration relation fo r our sample, compared to the 
model of lBullock et all ||2rj 



c{M,a) = K- 



(27) 



a c {M) 

where the collapse time a c is defined as the epoch at which 
the typical collapsing mass, M*(a c ), equals a fixed fraction 



M vlr [h-> Mj 

Figure 4. Relation between virial mass and NFW concentration 
parameter c = r v i r /r s . Solid squares are used for relaxed systems, 
empty squares for minor mer gers and stars for m ajor mergers. 
Dashed line shows the model o flBullock etldlfeOOll) . Dotted lines 
indicate the expected one-sigma scatter. 



F of t he halo mass at epoch a. According to lBullock et alJ 
teOOll) . we set K = 3 and F = 0.001 in order to account for 
the most massive haloes. The scatter around the relation is 
fitted by K = 2 and K = 4.7. 

Our sample of dark matter haloes is consistent with the 
expected trend, although our cluster-size haloes seem to be 
slightly more concentrated than the theoretical model. On 
the other hand, major mergers display systematically lower 
concentrations, since their collapse time is very close to the 
present. However, the value of a c is not well defined in a 
merging system. In particular, during the first stages of the 
merger the profile corresponds to an old object, perturbed 
by an approaching satellite. It is only after virialisation when 
the profile relaxes to its final form, and a c increases accord- 
ingly. 

3.2 Angular momentum 

In addition to the radial mass distribution, the kinetic struc- 
ture of numerical haloes can offer interesting insights into 
the formation of galaxies and galaxy clusters. In particu- 
lar, we are interested in the specific angular momentum of 
dark matter particles in order to set the eccentricity in the 
analytical model. 

We separate the velocity field into a random component 
(i.e. velocity dispersion) and ordered rotation (i.e. average 
j). The velocity dispersion of a collisionless CDM halo is 
related to the tota l mass distribution by virt ue of the Jeans 
equation (see e.g. iBinnev fc TremaineMl987l) . For a spher- 
ically symmetric system with isotropic velocity dispersion, 
neglecting infall, 



1 d(pa r 
P 



GM 



dr 



(28) 



where p is the local density, ay is the radial velocity dis- 
persion, and M is the mass enclosed within radius r. An 
approximate velocity dispersion profile could be derived by 
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Figure 5. Phase-space density profiles of our haloes (dots). Error 
bars show the rms scatter of individual profiles around the aver- 
age. The best-fit 'universal' profile 1-tlH is plotted as a dashed 
line. 



Figure 6. Contribution of bulk rotation (solid squares) and ran- 
dom motions (empty squares) to the tangential velocity of dark 
matter particles. Points represent the average over all haloes, and 
error bars indicate the rms scatter. 



substituting a given mass distribution (e.g. NFW) and set- 
ting an arbitrary normalisation oy(0). The contribution of 
random motions to the angular momentum of the CDM par- 
ticles would be given by the tangential velocity dispersion. 
In the isotropic case, this amounts to (j' 2 ') = 2r 2 a 2 (r). Note 
that random motions do not contribute to the total angular 
momentum of the halo (i.e. (J) = 0). 

A more empirical app roach to the velocity disp ersion 
profile has been followed bv lTavlor fc Navarro! i200fl) . They 
realised that the coarse-grained phase-space density of a 
sample of numerical galaxy-size haloes followed a power law 



(29) 



P 

^txr 
<j 

with /3 = —1.875 over more than two and a half decades 
in radius. Ras ia et alJ J2003) obtained a similar result for 
cluster-size haloes, although their best-fitting slope is f3 = 
— 1.95. The normalisation, though, has not been given in any 
of the two studies. 

We have investigated the phase-space structure of the 
CDM haloes in our sample, including both galaxies and 
galaxy clusters. The average profile p/cr 3 is plotted in Fig- 
ure^] We find that all our haloes are well described by 

-1.90±0.05 



= 10 



1.46±0.04 Pc 



(30) 



where i> vir = GM v i r /r v i T . T his result is in fair agree ment 
with the slopes re ported by iTavlor fc Navarrol fcOOll) and 
iRasia et al.l l)2003l) . The scatter around the average profile 
is remarkably low, taking into account that our haloes span 
four orders of magnitude in mass, and that we considered 
all systems (even major mergers) in the analysis. Indeed, 
we do not find any evidence that the slope of the relation 
depends on mass or dynamical state. Therefore, we claim 
that expression 1301 can be regarded as a 'universal' phase- 
space density profile with only one free parameter, which is 
the virial radius of the halo. 

We now show in Figure that angular momentum is 
indeed dominated by random motions. We compare the 



contribution to the tangential velocity of CDM particles 
from < j > and < j 2 > over spherical shells. The aver- 
age bulk rotation velocity of our halos (solid symbols in 
Figure is about O.ltVir, although we find a significant 
increase near the centre. Never the less, it is worth men- 
tioning that the separation in bulk and random velocity is 
prone to numerical errors at the innermost regions. The spe- 
cific angular momentum grows roughly linearly with radius 
(rigid bod y rotation), in agreement with previous numer- 
ical work jBarnes fc Efstathioul Il987t iBullock et alJ l200ll: 
Ivan den Bosch et^dTl2o"o2Tlciien'fe^ngll2'o02l) . In any case, 
the tangential velocity of individual dark matter particles is 
always dominated by the velocity dispersion, for all haloes, 
at any radius. 

As was discussed in Section 12.31 the orbit eccentricity 
is an important ingredient of our analytical model, since 
it measures how deep a CDM particle is able to penetrate 
within the gravitational potential. The ideal case e = cor- 
responds to a configuration in which all particles would orbit 
at constant radius with no mixing between concentric shells, 
while e = 1 (often assumed in spherical infall models) im- 
plies that every particle goes through the very centre of mass 
of the halo. For large values of the eccentricity, the central 
regions are mainly composed of particles coming from the 
outer shells, that happen to be near the pericentre of their 
orbits. 

The eccentricity profile of our haloes is plotted in 
Figure For each object, we compute the spherically- 
symmetric potential derived from its mass distribution. The 
pericentric and apocentric radii of every CDM particle are 
estimated from its position and velocity, and the eccentricity 
is computed as 

(31) 



' max ' mm 

The general trend is that particle orbits are slightly 
more radial as we move out to the current turn-around ra- 
dius of the halo, Rtn- Moreover, we find a systematic depen- 
dence on the dynamical state. As can be seen in Figure |H| 
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Figure 7. Orbit eccentricity of CDM particles as a function of 
apocentric radius. Dotted line marks e = 0.5; dashed line, expres- 
sion 1321 . 




10- 3 0.01 0.1 1 10 3 0.01 0.1 1 10- 3 0.01 0.1 1 

Figure 8. Same as Figure 171 separating relaxed haloes (left), 
minor (middle) and major (right) mergers. 



major mergers are well described by constant eccentricity up 
to the virial radius , while relaxed systems are more consis- 
tent with a power-law profile. Minor mergers are somewhat 
in the middle. The average profile can be fitted by a power 
law, but the slope is shallower than for relaxed systems. 
A least-square fit to the relaxed population yields 



e(r max )^0.8(^)' 



(32) 



for r max < 0.1-Rta- We note that our approximation to com- 
pute the eccentricity breaks down at larger radii, since the 
mass distribution (and thus, the gravitational potential) are 
by no means static beyond the virial radius. 



4 MODEL RESULTS 

We now attempt to reproduce the density profiles of our sim- 
ulated haloes with the model described in Section|5| In order 
to keep the number of free parameters to a minimum, we as- 
sumed a constant eccentricity, e ~ 0.5, which corresponds to 
the average over the whole radial range. This value implies 
that particles are able to sink into the dark matter potential 



1 For our dark matter haloes, 
of 0.2 - 0.3. 



i /Rtn is typically of the order 



Rf [Mpc/h] 

Figure 9. Best-fitting values of the peak height, u, and the 
smoothing scale, Rf . Black squares represent relaxed haloes, open 
squares are used for minor mergers, and stars for major merging 
systems. Shaded areas indicate yJx 2 /{dof) ^ 0.12. Dashed lines 
are drawn at constant r v ; r . 



up to one third of their maximum radius, or, equivalently, 
that all mass within radius r comes from shells whose turn- 
around radius was ^ 3r. The effect of a power law profile 
will be considered in Section f4.2l 



4.1 Constant eccentricity 

Once the eccentricity is set to e = 0.5, the model has two free 
parameters, which describe the primordial density peak: its 
height, v, and the smoothing scale, Rf . These parameters set 
the mass and formation time of the halo. For a given mass, 
a higher peak on a smaller scale corresponds to an earlier 
formation time and a steeper density profile near the centre. 

For each halo, we performed a \ 2 minimisation of the 
mass profile. Best-fitting values of Rf and v a re plotted in 
Figure^ as well as the areas where \/ X 2 l{dof) ^ 0.12. Most 
cluster-size haloes can be described as very high-i^ peaks, 
while galaxies seem to have collapsed around less extreme 
fluctuations. This means that clusters are expected to form 
earlier than galaxies, in the sense that the seeds of their 
dark matter haloes are already in place at a higher redshift. 

An exception to this rule are major mergers. While some 
of these systems still have an early collapse time, some oth- 
ers have collapsed approximately at the present epoch. The 
first class corresponds to objects that have not relaxed yet, 
and their core still corresponds to the old halo. As relax- 
ation completes and substructure is erased, the best-fitting 
parameters move along the lines of constant mass in the 
v — Rf diagram. The smoothing scale rises sharply, and v 
decreases accordingly. 

As can be seen in Figure [5] there is never the less a 
certain degeneracy between Rf and v, which prevents a re- 
liable determination of the formation time. The exact value 
of the best-fitting parameters may vary within the shaded 
area, depending on the details of the fit. For instance, we 
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Figure 10. Circular velocity profiles, v\ = GM/r. Crosses correspond to the simulation data, while dotted lines are the best-fitting 
NFW model and solid lines the best-fitting spherical collapse model. Relaxed haloes, minor and major mergers are plotted on the top, 
middle and bottom panels, respectively. 



tend to get higher peaks on smaller scales as we give more 
weight to the inner parts of the profile. 

This can be understood when we compare the numerical 
data with the results of our spherical collapse model. We 
chose to fit the radial range 0.1 < r/r v i r < 1 in order to test 



the quality of the extrapolations, both towards the centre 
and to large radii. Individual circular velocity profiles are 
shown in Figure Hoi together with the best fits provided by 
our model and NFW formula. 

Although the mass distribution is fairly well described 
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Figure 11. Accuracy of NFW (open squares) and spherical col- 
lapse (solid) models. Symbols correspond to the average over 
haloes classified as relaxed (left panel), minor (middle) or major 
mergers (right). Error bars indicate one-er scatter of individual 
profiles, and vertical dotted lines show the radial range used for 
the fit. 

in general terms, the central density is usually underesti- 
mated by both models. When we fit the innermost parts, we 
are biased towards more concentrated distributions. Indeed, 
if only data within O.lrvir are considered, the best-fitting 
profiles are rather unrealistic. 

At large radii, the spherical collapse model reproduces 
the upturn in the circular velocity, while NFW drops to 
zero. This is due to the fact that the density in the NFW 
model vanishes at infinity, while it tends asymptotically to 
the mean value in the spherical collapse model. Adding a 
constant density background to the NFW formula is enough 
to bring the profile in agreement with the numerical data. 

We assess the accuracy of both NFW and spherical col- 
lapse models in Figure ITT1 where the quantity 

AM = M mode l - Mdata . , 

M ~ M data [66) 

is plotted as a function of radius. 

The most important difference between both models is 
that the accuracy of NFW fit improves significantly for re- 
laxed systems, particularly concerning the extrapolation to- 
wards r — > 0. The fits based on our spherical collapse model 
are on average less accurate, and the scatter between in- 
dividual dark matter haloes is a little bit larger. Yet, the 
uncertainty is always lower than 20 per cent for relaxed ob- 
jects, and only slightly larger for merging systems, where 
the prediction of the spherical collapse model is indeed very 
similar to the NFW fit. 

4.2 Variable eccentricity 

Although constant eccentricity might provide a useful ap- 
proximation for merging systems, it is evident from Figure|H| 
that it is not a good description of relaxed haloes, where the 
eccentricity profile e(r) increases significantly from the cen- 
tre to the turn-around radius. Moreover, the assumption of 
constant eccentricity leads to systematic differences between 
the mass distribution predicted by our model and the nu- 
merical data, as shown in Figure [Til 

Therefore, we would like to investigate the consequences 
of using a more elaborate prescription for e(r). For constant 
eccentricity, high values of e lead to steeper density profiles 
in the central regions, whereas the outer parts of the halo 
remain largely unaffected (see Figure . We plot in Fig- 
ure 1121 the density profile resulting from our power- law fit 
{Hll for v = 3 and R f = 1 h' 1 Mpc. The mass distribution 
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Figure 12. Comparison between a NFW density profile (dot- 
ted line), the spherical collapse model with constant eccentricity 
e = 0.5 (dashed line), and a power-law e(r) given by 1321 (solid 
line). Top panel displays the density profile, and bottom panel its 
logarithmic slope. 
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Figure 13. Same as Figure ITTI when the eccentricity is set ac- 
cording to 1321 . 

is very similar to that obtained with e = 0.5, but the shape 
is slightly different. There is a small increase in the density 
between 10 and 100 kpc, but the profile flattens near 
the centre due to the more circular orbits. 

This flattening can be clearly seen on the bottom panel 
of Figure 1121 where we plot the logarithmic slope of the 
density profile. It is interesting to note that the spherical 
collapse model predicts a finite density at r — when non- 
radial motions are included, albeit the size of the 'core' is 
extremely small. As was discussed in Section 13.11 a similar 
trend is shown by the relaxed haloes in our sample. 

At our resolution limit, the net effect of using expres- 
sion 1321 is just to increase the density at small radii. Then 
our model gives a somewhat poorer description of merging 
systems, but the quality of the fit improves considerably for 
relaxed haloes. As can be seen in Figure FHfl the accuracy of 
our model is comparable to NFW formula when a realistic 
prescription is used for the eccentricity of particle orbits. 
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5 DISCUSSION AND CONCLUSIONS 

The radial density profile of dark matter haloes has been 
investigated within the framework of the spherical collapse 
theory. We have shown that the model described in Section[2] 
is able to reproduce the mass distribution of realistic CDM 
haloes. Although the final profile cannot be cast in a simple 
analytical form, it provides not only a mere phenomenologi- 
cal fit, but a physically-motivated description of the density 
distribution in terms of the primordial initial conditions. 

However, it is not easy to understand how the assump- 
tion of spherical symmetry could be able to describe the 
hierarchical assembly of cosmological structures. Instead of 
continuous infall of spherical shells, the formation of CDM 
haloes observed in numerical experiments takes place in a 
discrete and anisotropic way. Most of the matter is accreted 
in clumps, along the preferred directions set by the filamen- 
tary large scale structure. 

None the less, the very complica ted coalescence pro - 
cess looks very regular in ener gy space l|Zaroubi et all l996'l. 
Moreover. lMoore et alJ Jl999t) have shown that the final den- 
sity profile is not very sensitive to the details of the merging 
history, by comparing the mass distribution of a galaxy clus- 
ter halo with a re-simulation in which the power spectrum 
was truncated at ~ 4 Mpc scales. 

We suggest, as a plausible explanation for the success of 
the spherical infall model, that merging is implicitly taken 
into account through the smoothing scale R{. Cosmological 
structures do not form around maxima of primordial density 
field, but of the smoothed density field ( Hoff man fc Shahaml 
1985). Some memory of the initial conditions will be lost 
during major mergers, particularly at the innermost regions 
of the resulting halo. Contrary to the common view (see e.g. 
BBKS), we argue that Rf has a precise physical interpre- 
tation; below the mass scale defined by Rt, all information 
about the primordial substructure would have been erased 
by relaxation processes. 

In the outer regions, matter is accreted in a more gentle 
way. Minor mergers do not significantly alter the dynamical 
structure of the halo. The mass of infalling clumps is much 
less than the average over a spherical shell at large radii. 
Thus, the density profile and the accretion rate are deter- 
mined by the amount of matter available to the halo, which 
is ultimately set by the primordial initial conditions. 

Never the less, there is still an additional degree of free- 
dom, related to the amount and distribution of angular mo- 
mentum within the dark matter halo. Angular momentum 
sets the shape of the density profile at the inner regions. For 
pure radial orbits, the core is dominated by particles from 
the outer shells. As the angular momentum increases, these 
particles remain closer to their maximum radius, resulting 
in a shallower density profile. 

We found that angular momentum is dominated by the 
tangential component of the velocity dispersion. Random 
motions are well described by a 'universal' phase-space den- 
sity profile over several orders of magnitude in radius. This 
profile is a power law with slope ~ —1 .9, in agreement 
with the re sults oflTavlor <^ Nav arrol (|200lT l for galactic-size 
haloes and iRasia et £dTl2003ll for clusters. We found that 
this profile is valid for all objects in our sample (spanning 
a mass range from 10 11 to 10 15 M©) as long as the virial 
radius is used as a scale parameter. 



Although it would be desirable to establish a link be- 
tween the initial conditions and the specific angular mo- 
mentum, we resorted to a phenomenological approach. Ac- 
cording to the results of our numerical simulations, the or- 
bit eccentricity can be set to a constant value of e = 0.5 
as a first order approximation. However, relaxed haloes are 
better described by a power-law profile in terms of their 
turn-around radius. The details of the mass distribution de- 
pend only moderately on the prescription assumed for the 
angular momentum. Therefore, the spherical collapse model 
provides a valuable tool to predict the structure of virialised 
dark matter haloes, given the power spectrum of primordial 
fluctuations. Indeed, we claim that the physical origin of the 
mass distribution observed at the present day is related to 
the shape of the primordial density peaks. 'Universal' pro- 
files with two free parameters arise naturally from Gaussian 
random peak statistics, since the primordial fluctuations are 
fully specified by their height, i/, and smoothing scale, Rt. 
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